%% agro yield calculation (see 
q = param.q;
d = param.d(:,2);
qCS = cumsum(q);
dCS = cumsum(d);
figure(100);
plot(cumsum(q))
hold on,
plot(cumsum(d),'r') 
hold off
%%
fracQD = qCS./dCS;
figure(102); 
subplot(1,2,1), 
plot(cumsum(q))
hold on,
plot(cumsum(d),'r') 
hold off; xlabel('time'); ylabel('q and d');
legend('Avail. crop H_{2}O (q)','Tot. crop H_{2}O requirement d')
subplot(1,2,2);
semilogy(fracQD,'k'); hold on;
minFQD = min(fracQD);
iMinFQD = find(minFQD==fracQD);
plot(iMinFQD,minFQD,'ro'); xlabel('time'); ylabel('q/d');
legend('q/d','critical crop stage')
%% then
ky = .85;
yrc = zeros(param.years * param.nDays,1);
for t = 1 : param.nDays : (param.years-1) * param.nDays
% for t = 1 : param.years * param.nDays
    yrc(t:t+param.nDays) = min(1 - ky * (1 - cumsum(q(t:t+param.nDays))' ./ cumsum(d(t:t+param.nDays))), ...
        1 - ky * (1 - sum(q(t:t+param.nDays)) / sum(d(t:t+param.nDays))));
    minYRC(t) = min(yrc(t:t+param.nDays));
    iMinYRC(t) = find(minYRC(t)==yrc(t:t+param.nDays),1)+t;
end
figure(101);
plot(yrc,'k'); 
hold on;
scatter(iMinYRC,minYRC,'ro');
% plot(1-yrc), this then is going to be the reward for RL
hold off
%% yield derived from end of growing season harvest
ym = zeros(size(q));
ym(200:250) = 1;
ym = ym .* minYRC .* d';
plot(cumsum(ym)) 